This is a pre-submission DRAFT. Some of this text will be moved into the Methods and Results section of the manuscript’s main text. Other text may remain in supplemental materials. Details that may be important now, but should likely be removed before publication are given as block quotes. While any process is running on nscc, you can monitor it’s progress indirectly by checking the CPU usage on NSCC node 28.

In order to characterize the development of polyphenic traits in the wings, flight muscles and gonads of the soapberry bug, Jadera haematoloma, we sequenced transcriptomes from a multidimensional matrix of life stage, tissue, sex, food regime and host population.

Sampling design

Transcriptome studies proceeded through three phases. Initially, we sampled whole bodies from 12 individuals, including three biological replicates of each sex and morph combination from the Plantation Key, FL population (phase 1). Next we sampled dorsal thorax and gonads from three individuals of each sex and morph combination from the Plantation Key, FL and Aurora, CO populations (phase 2). Samples from these phases were prepared as full-length 150-bp paired-end Illumina libraries and used to assembly a reference transcriptome.

The final phase utilized 3’-tag sequencing with the exclusive goal of quantifying gene expression. The sampling design at this phase included two stages (fifth instars and nascent adults), two tissues (dorsal thorax and gonads), two sexes, long- and short-wing adult morphs, high and low food regimes, four populations (collected from Key Largo, FL, Plantation Key, FL, Aurora, CO, and Frederick, MD), and three biological replicates of each unique combination of the preceding factors. Thus the full matrix contained 288 samples. Seven samples were lost or did not pass library quality controls, resulting in 281 samples in the analysis.

Shape analysis

Most J. haematoloma can be qualitatively identified as “long winged” or “short winged” based on the shape of the wings and whether their distal tips extend beyond the posterior of the abdomen when folded up at rest. However, closer inspection reveals that some individuals have intermediate shapes and lengths. We reasoned that gene expression was likely to also follow a more continuous distribution. Therefore, rather than relying on a categorical assignment of specimens into morphs, gene expression was modeled based on the quantification of wing shapes. (Categorical modeling by morph was also performed, and the two approaches identified similar lists of top DEGs.) Thoracic muscles are required for flight, therefore we also tested for correlations between gene expression and thorax shape.

Bugs were anesthetized using CO2 and imaged on a trinocular stereo microscope (VWR International, Radnor, Pennsylvania, USA) with a Moticam 5 digital camera. Dorsal and ventral images of each specimen were recorded with a millimeter-scale ruler to scale pixel measurements. ImageJ was used to place landmarks.

Adult wing shape was quantified as described by Fawcett et al. (2018).

Table SX: Wing landmark positions

landmark position
1 Bifurcation of the radial and medial veins
2* Maximum curvature of clavus, near the posterior tip of the mesonotum
3-4* Two equally spaced semilandmarks along the claval margin, anterior of landmark 2
5 Bifurcation of the radial vein and R2 crossvein
6 Intersection of the R2 crossvein and the medial vein
7 Intersection of the cubitus vein and the first anal vein
8* One semilandmark along the claval margin, between landmarks 2 and 7
9 Intersection of the medial vein and the posterior crossvein
10 Intersection of the radial vein and the anterior crossvein
11 Terminus of the radial-medial crossvein
12 First branch of the anterior membrane vein
13* Maximum posterior-distal curvature of the membrane region of the wing (roughly lateral to landmark 11 when held in place on the bug)
14-16* Three equally spaced semilandmarks along the posterior wing margin, between landmarks 7 and 13
17-21* Five equally spaced semilandmarks along the distal wing margin, between landmarks 11 and 13
22-24* Three equally spaced semilandmarks along the anterior wing margin, between landmarks 1 and 11
* semilandmark

This system of landmarks was adapted to quantify juvenile wing pad shapes as described in Table SX.

Table SX: Juvenile wing pad landmark positions.

landmark position
1 Bifurcation of the radial and medial veins
2 Medial base of the wing pad
3-4* Two equally spaced semilandmarks along the base of the wing pad, at its attachment to the body, between landmarks 1 and 2
5 Bifurcation of the radial vein and R2 crossvein
6 Intersection of the R2 crossvein and the medial vein
7 Intersection of the cubitus vein and the first anal vein
8* One semilandmark along the claval margin, between landmarks 2 and 7
9 Intersection of the medial vein and the posterior crossvein
10 Intersection of the radial vein and the anterior crossvein
11 Terminus of the radial-medial crossvein
12* Maximum posterior-distal curvature of the wing (roughly lateral to landmark 11)
13* One equally spaced semilandmark along the posterior wing margin, between landmarks 7 and 12
14-18* Five equally spaced semilandmarks along the distal wing margin, between landmarks 12 and 11
19-21* Three equally spaced semilandmarks along the anterior wing margin, between landmarks 1 and 11
* semilandmark

Thorax shapes of juvenile and adult bugs were digitized using the same set of landmarks (Table SX) .

Table SX: Thorax landmark positions.

landmark position
1 Anterior right of the pronotum, as viewed from above
2 Posterior right of the pronotum
3 Posterior left of the pronotum
4 Anterior left of the pronotum
5 Anterior midline of the pronotum
6 Anterior left of the mesonotum
7 Anterior midline of the mesonotum
8 Anterior right of the mesonotum
9 Posterior of the mesonotum

Variation in shapes was analyzed using the geometric morphometric methods implemented in the R packages geomorph and borealis. Generalized Procrustes analysis with partial Procrustes superimposition was preformed using minimized bending energy. The shape of individual specimens was quantified along the first principle component axis (Fig. SX,SX). Procrustes ANOVA with permutation was used to assess hypotheses for patterns of shape variation among the aligned specimens, using 10,000 iterations.

Figure SX. Morphospace plotting the first two principle component axes for the shapes of juvenile wing pads (A) and adult wings (B). Specimens are color-coded by sex and morph (for adults), and colored convex hulls show the space occupied by these groups. Deformation grids illustrate shapes at the extremes of each axis. (Anterior is up; distal is right.) The percent of total shape variance described by each principle component in given on the axis.

Figure SX. Morphospace plotting the first two principle component axes for the shapes of juvenile (A) and adult pronotum and mesonotum (B). Specimens are color-coded by sex and morph (for adults), and colored convex hulls show the space occupied by these groups. Deformation grids illustrate shapes at the extremes of each axis. (Anterior is up.) The percent of total shape variance described by each principle component in given on the axis.

Nucleic acid extractions

Individuals were collected from laboratory cultures, anesthetized using CO2 exposure, photographed, and flash frozen in liquid nitrogen. Samples were stored at -80˚C until further processing. The thorax was removed from head and abdomen with a clean scalpel blade, and the dorsal thorax was separated above the legs while tissue was still frozen. Gonads were then dissected from the abdomen of adults. For juveniles, the entire abdomen was included. Nucleic acid extractions used the Invitrogen PureLink RNA extraction kit.

Sequencing and assembly of a reference transcriptome

Samples used for transcriptome assembly (phases 1 and 2) were delivered on dry ice to Beckman Coulter Genomics (Danvers, Massachusetts) for poly-A selection, preparation of 125 bp paired-end TruSeq libraries, and sequencing using Illumina HiSeq. Libraries from phases 1 and 2 were sequenced separately. Reads from all samples were trimmed using Trimmomatic version 0.33 (Bolger et al. 2014) with the parameters ILLUMINACLIP:/export/local/src/Trimmomatic-0.33/adapters/TruSeq3-SE.fa:2:30:10 HEADCROP:12 LEADING:3 TRAILING:3 SLIDINGWINDOW:3:15 MINLEN:50. The resulting reads were assembled into transcripts by Trinity version 2.0.6 (Grabherr et al. 2011). Before annotation, potentially redundant transcripts with greater than 90% identity were collapsed using CD-HIT-EST version 4.6 (Li & Godzik 2006).

cd /research/drangeli/phase2_BCG_RNAseq/Jhae_GU_trinity_assembly/
cd-hit-est -i GU.Trinity.fa -o GU.Trinity.c90.fa -c 0.90 -n 10 -d 0 -M 64000 -T 12

This took about 3.5 days on node 26.

TransDecoder identified coding sequences. At each stage, BUSCO version 5.2.2 (Manni et al. 2021) was used to assess transcriptome completeness.

busco --in GU.Trinity.fa \
--lineage_dataset /research/drangeli/DB/busco.lineages/hemiptera_odb10 \
--out BUSCO.report.for.GU \
--augustus --augustus_species fly --mode transcriptome --cpu 24

Consolidation by cd-hist-est dramatically reduced the number of sequences flagged as duplicates (Table S1).

Table S1: Details of transcriptome datasets at each stage of processing.

dataset contigs BUSCO short summary
initial assembly 549,485 C:97.0%[S:54.2%,D:42.8%],F:1.3%,M:1.7%,n:2510
initial assembly CDS 85,937
90% consolidation 395,125 C:96.6%[S:77.7%,D:18.9%],F:1.4%,M:2.0%,n:2510
90% consolidation CDS 50,771 C:94.4%[S:77.8%,D:16.6%],F:2.2%,M:3.4%,n:2510

Annotation of the transcriptome

Transcript annotation was performed by EnTAP (Hart et al. 2019) using the NCBI invertebrate RefSeq dataset, SwissProt and Trembl databases. EnTAP assigns sequence names and GO terms based on similarity.

EnTAP --runP --threads 24 \
-i /research/drangeli/phase2_BCG_RNAseq/Jhae_GU_trinity_assembly/GU.Trinity.c90.fa \
-d /export/groups/drangeli/entap_config/bin/invert.dmnd \
-d /export/groups/drangeli/entap_config/bin/uniprot_sprot.dmnd \
-d /export/groups/drangeli/entap_config/bin/uniprot_trembl.dmnd \
--ini /export/groups/drangeli/Jhae.genome/entap_config.ini

The J. haematoloma sequences were matched to other Hemiptera, primarily Halyomorpha halys, Cimex lectularius and Nilaparavata lugens for the RefSeq comparisons and Oncopeltus, Lygus and Riptortus from TrEMBL. The contaminant sequence species are also as expected, coming from bacteria and yeasts. Only 12 potential contaminant sequences had homology to humans.

Sequencing for gene expression quantification

For gene expression analysis, RNA samples were shipped overnight on dry ice to the DNA Technologies Core at the University of California at Davis. Libraries were prepared using the QuantSeq 3’ mRNA-Seq Library Prep Kit (Lexogen, Greenland, New Hampshire, USA) in the three batches for sequencing in separate Illumina HiSeq lanes. Batch 3 was sequenced twice and reads counts for each run were added together for each sample after filtering.

These batches were sequenced in three lanes yielding 115.9×109 bp. These data

Short read sequences were inspected for quality using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), before and after trimmed using Trimmomatic version 0.33, as in the example command below.

java -jar trimmomatic-0.33.jar SE \
raw.data/A1-d25-2_S161_L006_R1_001.fastq \
filtered.reads/A1-d25-2_S161_L006_R1_001.trimmed.fastq \
ILLUMINACLIP:/export/local/src/Trimmomatic-0.33/adapters/TruSeq3-SE.fa:2:30:10 \
HEADCROP:16 LEADING:3 TRAILING:3 SLIDINGWINDOW:3:15 MINLEN:79

A simple bash script was used to apply trimming to all read files and save the output in a separate folder.FastQC inspection after trimming revealed notable improvements.

Table S2: Read numbers before and after quality filtering.

batch raw reads filtered reads passing filter
1 341,370,316 323,188,378 94.7%
2 321,564,969 313,590,249 97.5%
3 484,588,962 419,931,564 86.7%
overall 1,147,524,247 1,056,710,191 92.1%

Mapping 3seq reads to the transcriptome

Kallisto was used to quantify the abundance of transcripts from 3seq data using a pseudo-alignment method (Bray et al. 2016). First, a Kallisto index was created from the assembled transcripts, collapsed by cd-hit-est.

We chose to use this target, rather than the TransDecoder CDS files, because if there are fragmentary or non-coding RNAs that were sequenced, then some reads will correspond to those RNAs too.

cd /research/drangeli/phase2_BCG_RNAseq/Jhae_GU_trinity_assembly/
kallisto index --index=GU.Trinity.c90.kallisto.index GU.Trinity.c90.fa

In order to run Kallisto on single-end sequence, the program requires the mean and standard deviation of read lengths. The total number of reads was found using bash commands.

cd /export/groups/drangeli/morphDE/davis.3seq/filtered.reads

gunzip -c 180824/*.fastq.gz 190917/*.fastq.gz 200611/*.fastq.gz 200701/*.fastq.gz \
| sed -n -e '2~4p' \
| awk '{ print length($0); }' > read.lengths.txt

wc -l read.lengths.txt

Total read number was 1,149,235,725.

Mean and standard deviation were calculated using awk.

cat read.lengths.txt | awk '{for(i=1;i<=NF;i++) {sum[i] += $i; sumsq[i] += ($i)^2}} END {for (i=1;i<=NF;i++) { printf "%f %f \n", sum[i]/NR, sqrt((sumsq[i]-sum[i]^2/NR)/NR)} }'

A bash script was used to run kallisto quant for each sample and to rename the output files with the sample names.

The result is a folder kallisto.vs.JhaeGUtx that contains .abundance.tsv files for each sample, listing the counts for each transcript.

The resulting output files were then merged into a single comma-separated values (CSV) file with columns for each sample.

This file kallisto.counts.csv is 309 Mb, which is too large to post to GitHub.

Differential gene expression analysis

Gene expression differences were examined using the Bioconductor package DESeq2 version 1.32.0 (Love et al. 2014) in R version 4.1.1 (R Core Team 2021). The function DESeq fits a generalized linear model where counts are modeled using a negative binomial distribution with a gene-specific dispersion parameter. Means were estimated using approximate posterior estimation for generalized linear models (Zhu et al. 2018). We extracted log2 fold-change and p-values corrected by independent hypothesis weighing (Ignatiadis et al. 2016) and applied a critical threshold of 0.05. Contrasts were made for sex, wing morph and food regime, after accounting for batch effects (multiple sequencing runs). Data were analyzed separately for each stage (fifth instar juveniles and nascent adults) and for each tissue (thorax and gonads). Differential gene expression was visualized by customized plots generated using ggplot2. Variance stabilizing transformation, using the function vst, was applied before ordination using principle component analysis.

For more details see the R script analysis.dge.R

Batch effects were introduced by the need to sequence samples in multiple lanes (Figure S1).

Figure S1. Principle components of uncorrected counts, color-coded by biological factors (A) and sequencing batch (B).

DESeq2 allows variance from confounding factors, such as sequencing batch, to be accounted for through generalized linear modeling. The normalized counts that result from the model then reflect biologically meaningful variance (Figure S2).

Figure S2. Principle components of transcript counts for each sample, following normalization for sequencing run (batch). Samples are color-coded by biological origin.

More complex models can be applied to examine specific contrasts of biological interest. These models partition the variance in gene expression sequentially to factors in the model. For more complex models, less power exists to detect differences, but it allows for the influence of multiple factors to be disentangle from one another.

Expression differences by sex

Sex is often an important factor effecting gene expression. While it is not our focus of study, examining our dataset for differences between females and males in each tissue is a useful confirmation of these expectations (Figure S3. Moreover, sexually dimorphic gene expression has not previously been described in J. haematoloma. Genes were differentially expressed at both sampled stages and in both sampled tissues. Adults had a larger number of sexually dimorphic DEGs (84% in gonad; 29% in thorax) compared to juveniles (6.1% in gonads, 0.18% in thorax). The gonads also had more sexually dimorphic gene expression than did the thorax. Within the gonads, males had more significantly biased genes than did females, at both stages (Figure S3C-D). Several genes with roles in somatic sex determination were among these, including doublesex (dsx) and transformer [ref]. The histone acetyltransferase encoded by mof is known to increase expression on the X chromosome in D. melanogaster males (Smith et al. 2000), although the mof ortholog appears up-regulated in the J. haematoloma female thorax (Fig S3B). In the adult gonad contrast, upd is significantly biased towards male expression. In D. melanogaster testis, upd is expressed by hub cells promote maintenance of spermatogenic stem cells (reviewed by Matunis et al 2012).

Figure S3. Volcano plots, showing log2 fold-change (LFC) and -log10 FDR-adjusted p-values for contrasts by sex in the juvenile thorax (A), adult thorax (B), juvenile gonads (C) and adult gonads (D). Genes are highlighted in red if they have an absolute LFC values > 1 (two-fold difference) and an FDR-adjusted p-value < 0.05 (panels A) or 10-3 (panels B-D). Highlighted transcript are annotated where possible.

Expression differences by morph

Add here!

Figure S4. Volcano plots, showing log2 fold-change (LFC) and -log10 FDR-adjusted p-values for contrasts by morph in the adult thorax (A), adult gonads (B), adult ovaries (C) and adult testes (D). Genes are highlighted in red if they have an absolute LFC values > 1 (two-fold difference) and an FDR-adjusted p-value < 10-3 (panels A) or 0.05 (panels B-D). Highlighted transcript are annotated where possible.

Expression differences by stage

Add here!

Figure S5. Volcano plots, showing log2 fold-change (LFC) and -log10 FDR-adjusted p-values for contrasts by stage in the thorax (A), ovaries (B), and testes (C). Genes are highlighted in red if they have an absolute LFC values > 1 (two-fold difference) and an FDR-adjusted p-value < 10-3. Highlighted transcript are annotated where possible.

Expression differences by food regime

Add here!

Figure S5. Volcano plots, showing log2 fold-change (LFC) and -log10 FDR-adjusted p-values for contrasts by food regime in the juvenile thorax (A), adult thorax (B), juvenile gonads (C), and adult gonads (D). Genes are highlighted in red if they have an absolute LFC values > 1 (two-fold difference) and an FDR-adjusted p-value < 0.05. Highlighted transcript are annotated where possible.

Enrichment analysis

Enrichment analyses were performed based on annotations with terms from the GO and KEGG ontologies, as well as EggNOG protein domains, applied to transcripts by EnTAP. We also tested for enrichment of xenic or contaminant transcripts among those that were differentially expressed, and neither of these categories was significantly over-represented. For each term, p-values from the hypergeometric test, implemented in the base R function phyper, were adjusted for false discovery rates (FDR).

For more details on the enrichment analysis see the R script enrichment.R

Table S3: Terms from GO and KEGG ontologies and EggNOG protein domains over-represented in differentially expressed genes in each dataset and contrast.

stage tissue model factors ontology enriched terms top terms
juvenile thorax sex GO 0
juvenile thorax sex KEGG 0
juvenile thorax sex eggNOG 15 dsx/mab-3 DNA-binding domain, 9.13x10-7; Complementary sex determiner-like SR domain, 0.0180; proline-rich protein BAT2 N-terminal domain, 0.0180; CD34/Podocalyxin family, 0.0180; DFDF domain, 0.0180; Sex determination protein N terminal domain, 0.0180; Barentsz/CASC3 domain, 0.0215; Meckelin, 0.0215; CDC48, 0.0298
adult thorax sex GO 267 GO:0022626(L=5) cytosolic ribosome, 2.98×10^-41; GO:0006412(L=5) translation, 1.07×10-29; GO:0003723(L=4) RNA binding, 1.08×10-27; GO:0034645(L=4) cellular macromolecule biosynthetic process, 4.73×10-27; GO:0009059(L=4) macromolecule biosynthetic process, 1.16×10-26; GO:0000022(L=5) mitotic spindle elongation, 4.66×10-26; GO:0051231(L=4) spindle elongation, 3.81×10-25; GO:0044391(L=5) ribosomal subunit, 5.19×10-25; GO:0022625(L=6) cytosolic large ribosomal subunit, 9.80×10-24; GO:0044249(L=3) cellular biosynthetic process, 3.9×10-23
adult thorax sex KEGG 3 map04974 Protein digestion and absorption, 6.71x10-6; map04972 Pancreatic secretion, 1.04x10-4; map00513 Various types of N-glycan biosynthesis,7.60x10-3
adult thorax sex eggNOG 26 Reverse transcriptase, 2.11x10-9; Endonuclease/Exonuclease/phosphatase family, 1.38x10-8; Rnase H, 7.70x10-7; ZnF C2HC, 5.78x10-6; Lipoprotein N-terminal domain, 6.48x10-5; Vitellogenin-like N-terminal domain, 6.48x10-5; DUF1943, 2.54x10-4; Co/Zn superoxide dismutase, 8.63x10-4; CUB domain, 1.77x10-3; Trypsin, 3.17x10^-3
juvenile gonad sex GO 1299 GO:0005622(L=3) intracellular anatomical structure, 1.40×10-94; GO:0043229(L=3) intracellular organelle, 7.32×10-84; GO:0043231(L=4) intracellular membrane-bounded organelle, 1.47×10-64; GO:0043227(L=2) membrane-bounded organelle, 5.40×10-64; GO:0005737(L=4) cytoplasm, 2.40×10-55; GO:0044237(L=2) cellular metabolic process, 7.05×10-55; GO:0044238(L=2) primary metabolic process, 1.08×10-52; GO:0016043(L=2) cellular component organization, 3.73×1-51; GO:0071704(L=2) organic substance metabolic process, 6.15×10-50; GO:0043228(L=2) non-membrane-bounded organelle, 1.66×10-44
juvenile gonad sex KEGG 31 map01100 Metabolic pathways, 4.58×10-9; map05012 Parkinson disease, 1.22×10-7; map00710 Carbon fixation in photosynthetic organisms, 1.83×10-6; map05010 Alzheimer disease, 3.77×10-6; map00051 Fructose and mannose metabolism, 9.74×10-5; map01110 Biosynthesis of secondary metabolites, 1.04×10-4; map00030 Pentose phosphate pathway, 9.86×10-4; map04340 Hedgehog signaling pathway, 9.86×10-4; map05152 Tuberculosis, 9.86×10-4; map05110 Vibrio cholerae infection, 1.13×10-3
juvenile gonad sex eggNOG 46 COIL, 1.62x10-69; S/T kinase, 5.50x10-9; TRANS, 3.37x10-8; WD40 domain, 3.31x10-7; RNA recognition motif, 5.60x10-7; protein kinase, 1.14x10-6; Protein tyrosine kinase, 1.14x10-6; RNA recognition motif, 1.28x10-6; M17 cytosol aminopeptidase, 2.15x10-5; Leucin-rich repeats, 4.30x10-5
adult gonad sex GO 0
adult gonad sex KEGG 0
adult gonad sex eggNOG 0
juvenile thorax sex + wing pad shape GO NA (no DEGs)
juvenile thorax sex + wing pad shape KEGG NA (no DEGs)
juvenile thorax sex + wing pad shape eggNOG NA (no DEGs)
adult thorax sex + wing shape GO 1137 GO:0005737(L=4) cytoplasm, 1.73×10^-129; GO:0005622(L=3) intracellular anatomical structure, 8.29×10^-113; GO:0005739(L=5) mitochondrion, 7.11×10^-105; GO:0043229(L=3) intracellular organelle, 1.04×10^-93; GO:0043227(L=2) membrane-bounded organelle, 1.66×10^-79; GO:0043231(L=4) intracellular membrane-bounded organelle, 2.86×10^-79; GO:0044237(L=2) cellular metabolic process, 4.08×10^-67; GO:0005759(L=6) mitochondrial matrix, 1.58×10^-55; GO:0071704(L=2) organic substance metabolic process, 7.32×10^-55; GO:0044238(L=2) primary metabolic process, 1.36×10^-53
adult thorax sex + wing shape KEGG 15 map01100 Metabolic pathways, 2.50E-33; map05012 Parkinson disease, 4.47E-28; map05010 Alzheimer disease, 1.91E-19; map01110 Biosynthesis of secondary metabolites, 2.05E-12; map04974 Protein digestion and absorption, 4.23E-04; map01120 Microbial metabolism in diverse environments, 1.54E-03; map00640 Propanoate metabolism, 1.54E-03; map04260 Cardiac muscle contraction, 4.27E-03; map04972 Pancreatic secretion, 6.61E-03; map00513 Various types of N-glycan biosynthesis, 8.26E-03
adult thorax sex + wing shape eggNOG 54 COIL, 3.39E-41; TRANS, 1.20E-38; SIGNAL, 8.55E-17; Immunoglobulin I-set domain, 1.95E-04; Immunoglobulin C-2 Type, 2.37E-04; EF-hand, calcium binding motif, 6.37E-04; Immunoglobulin V-set domain, 7.20E-04; Tubulin, 9.85E-04; CUB domain, 1.23E-03; Immunoglobulin domain, 1.57E-03
juvenile gonad sex + wing pad shape GO 0
juvenile gonad sex + wing pad shape KEGG 3 map04975 Fat digestion and absorption, 0.00460; map00561 Glycerolipid metabolism, 0.0106; map04972 Pancreatic secretion, 0.0180
juvenile gonad sex + wing pad shape eggNOG 1 Lipase, 0.00992
adult gonad sex + wing shape GO 99 GO:0004197(L=6) cysteine-type endopeptidase activity, 0.0007786; GO:0045169(L=5) fusome, 0.0007786; GO:0000323(L=6) lytic vacuole, 0.0008834; GO:0005764(L=7) lysosome, 0.0008834; GO:0008219(L=3) cell death, 0.0008834; GO:0035070(L=5) salivary gland histolysis, 0.0008834; GO:0035071(L=6) salivary gland cell autophagic cell death, 0.0008834; GO:0048102(L=5) autophagic cell death, 0.0008834; GO:0008234(L=5) cysteine-type peptidase activity, 0.001606; GO:0005773(L=5) vacuole, 0.001743
adult gonad sex + wing shape KEGG 5 map04612 Antigen processing and presentation, 4.35e-6; map04145 Phagosome, 1.19e-5; map05205 Proteoglycans in cancer, 2.47e-5; map05323 Rheumatoid arthritis, 3.33e-3; map04142 Lysosome, 1.72-2
adult gonad sex + wing shape eggNOG 12 Peptidase C1, 1.28E-06; Cathepsin propeptide inhibitor domain (I29), 1.67E-06; SIGNAL, 0.00337; Thiol-activated cytolysin, 0.00952; Lipase 2, 0.0130; TRANS, 0.0157; Transmembrane amino acid transporter, 0.0253; Lipase, 0.0262
adult ovaries sex + wing shape GO 31 GO:0010447(L=4) response to acidic pH, 0.01403; GO:0030682(L=9) mitigation of host defenses by symbiont, 0.01403; GO:0052173(L=5) response to defenses of other organism, 0.01403; GO:0052200(L=6) response to host defenses, 0.01403; GO:0052572(L=7) response to host immune response, 0.01403; GO:0071291(L=6) cellular response to selenium ion, 0.01403; GO:0075136(L=5) response to host, 0.01403; GO:0009268(L=3) response to pH, 0.01904; GO:0010269(L=4) response to selenium ion, 0.01904; GO:0008121(L=5) ubiquinol-cytochrome-c reductase activity, 0.03328
adult ovaries sex + wing shape KEGG 0
adult ovaries sex + wing shape eggNOG 6 MIG-14 Wnt-binding factor, 7.11E-03; Fe/Mn superoxide dismutase, 7.11E-03; Pacifastin serine protease inhibitor, 0.0186; Reverse transcriptase, 0.0466; Pao retrotransposon peptidase, 0.0476
adult testes sex + wing shape GO 156 GO:0000131(L=4) incipient cellular bud site, 0.02493; GO:0000148(L=5) 1,3-beta-D-glucan synthase complex, 0.02493; GO:0000935(L=4) division septum, 0.02493; GO:0001837(L=6) epithelial to mesenchymal transition, 0.02493; GO:0004028(L=5) 3-chloroallyl aldehyde dehydrogenase activity, 0.02493; GO:0005773(L=5) vacuole, 0.02493; GO:0005934(L=3) cellular bud tip, 0.02493; GO:0005937(L=4) mating projection, 0.02493; GO:0006068(L=7) ethanol catabolic process, 0.02493; GO:0007431(L=5) salivary gland development, 0.02493
adult testes sex + wing shape KEGG 0
adult testes sex + wing shape eggNOG 21 Acyl-CoA reductase (LuxC), 0.0185; lysis protein, 0.0185; Alpha-lactalbumin / lysozyme C, 0.0185; asparaginyl peptidase C13, 0.0185; SIGNAL, 0.0185; Transglycosylase SLT domain, 0.0185; Thyroglobulin type-1 repeat, 0.0185; Transglycosylase-like domain, 0.0185; Thyroglobulin type I repeats, 0.0185; Cystatin-like domain, 0.0189
juvenile thorax sex + thorax shape GO NA (no DEGs)
juvenile thorax sex + thorax shape KEGG NA (no DEGs)
juvenile thorax sex + thorax shape eggNOG NA (no DEGs)
adult thorax sex + thorax shape GO 0
adult thorax sex + thorax shape KEGG 0
adult thorax sex + thorax shape eggNOG 22 SIGNAL, 3.13E-08; Calponin homology (CH) domain, 3.37E-04; TRANS, 9.67E-04; juvenile hormone binding protein (JHBP), 0.00285; Calponin, 0.00517; calmodulin-regulated spectrin-associated protein, 0.00517; CRAL/TRIO domain, 0.00517; Polysaccharide deacetylase, 0.00517; Sec14p-like lipid-binding domain, 0.00517; CRAL/TRIO N-terminal domain, 0.0118
juvenile gonad sex + thorax shape GO NA (no DEGs)
juvenile gonad sex + thorax shape KEGG NA (no DEGs)
juvenile gonad sex + thorax shape eggNOG NA (no DEGs)
adult gonad sex + thorax shape GO NA (no DEGs)
adult gonad sex + thorax shape KEGG NA (no DEGs)
adult gonad sex + thorax shape eggNOG NA (no DEGs)
juvenile thorax sex + food density GO 0
juvenile thorax sex + food density KEGG 0
juvenile thorax sex + food density eggNOG 15 DUF3480, 4.27E-03; Rad50 zinc hook motif, 4.27E-03; Smad anchor for receptor activation, 4.27E-03; DNA Topoisomerase I, 4.27E-03; Lipocalin, 7.11E-03; Triabin, 7.11E-03; von Willebrand factor type C domain, 0.0204; COIL, 0.0214; FYVE-like zinc finger, 0.0255; DnaJ/Hsp40, 0.0389
adult thorax sex + morph + food density GO 3 GO:0006801(L=4) superoxide metabolic process, 0.01508; GO:0072593(L=3) reactive oxygen species metabolic process, 0.03662; GO:0005319(L=3) lipid transporter activity, 0.04003; GO:0006869(L=4) lipid transport, 0.05938; GO:0010876(L=3) lipid localization, 0.05938; GO:0030710(L=4) regulation of border follicle cell delamination, 0.05938; GO:0033036(L=2) macromolecule localization, 0.05938; GO:0046956(L=6) positive phototaxis, 0.05938; GO:0047396(L=3) glycosylphosphatidylinositol diacylglycerol-lyase activity, 0.05938; GO:0052128(L=5) positive energy taxis, 0.05938
adult thorax sex + morph + food density KEGG 0
adult thorax sex + morph + food density eggNOG 15 DUF1943, 0.000260; Lipoprotein N-terminal Domain, 0.000260; Co/Zn superoxide dismutase, 0.000260; Vitellogenin, 0.000260; von Willebrand factor type D domain, 0.000389; L27 domain, 0.00579; SIGNAL, 0.00579; Chromatin associated protein KTI12, 0.00883; 6-phosphofructo-2-kinase, 0.00980; Guanylate kinase, 0.0138
juvenile gonad sex + food density GO 30 GO:0004044(L=5) amidophosphoribosyltransferase activity, 0.01636; GO:0004360(L=6) glutamine-fructose-6-phosphate transaminase (isomerizing) activity, 0.01636; GO:0007278(L=4) pole cell fate determination, 0.01636; GO:0016628(L=4) oxidoreductase activity, acting on the CH-CH group of donors, NAD or NADP as acceptor, 0.01636; GO:0019090(L=5) mitochondrial rRNA export from mitochondrion, 0.01636; GO:0019093(L=4) mitochondrial RNA localization, 0.01636; GO:0030719(L=5) P granule organization, 0.01636; GO:0030870(L=6) Mre11 complex, 0.01636; GO:0032440(L=5) 2-alkenal reductase [NAD(P)+] activity, 0.01636; GO:0070548(L=5) L-glutamine aminotransferase activity, 0.01636
juvenile gonad sex + food density KEGG 10 map04010 MAPK signaling pathway, 0.0195; map04141 Protein processing in endoplasmic reticulum, 0.0195; map04144 Endocytosis, 0.0195; map04612 Antigen processing and presentation, 0.0195; map05134 Legionellosis, 0.0195; map05145 Toxoplasmosis, 0.0195; map05162 Measles, 0.0195; map00520 Amino sugar and nucleotide sugar metabolism, 0.0195; map05164 Influenza A, 0.0226; map04915 Estrogen signaling, 0.0239
juvenile gonad sex + food density eggNOG 15 Rad50 zinc hook motif, 3.22E-03; Sugar isomerase domain, 3.22E-03; dwarfin domain B, 3.68E-03; Glutamine amidotransferases class-II, 3.68E-03; Interferon-regulatory factor 3, 3.68E-03; C-terminal Mad Homology 2 domain, 3.68E-03; Survival motor neuron protein, 3.68E-03; Collagen, 4.82E-03; Hsp70, 4.82E-03; MreB/Mbl domain, 4.82E-03
adult gonad sex + morph + food density GO 4 GO:0004446(L=7) inositol-hexakisphosphate phosphatase activity, 0.01454; GO:0052826(L=8) inositol hexakisphosphate 2-phosphatase activity, 0.01454; GO:0003993(L=6) acid phosphatase activity, 0.03929; GO:0052745(L=6) inositol phosphate phosphatase activity, 0.03929; GO:0006590(L=4) thyroid hormone generation, 0.06149; GO:0006768(L=4) biotin metabolic process, 0.06149; GO:0006797(L=4) polyphosphate metabolic process, 0.06149; GO:0007450(L=7) dorsal/ventral pattern formation, imaginal disc, 0.06149; GO:0008523(L=4) sodium-dependent multivitamin transmembrane transporter activity, 0.06149; GO:0015111(L=7) iodide transmembrane transporter activity, 0.06149
adult gonad sex + morph + food density KEGG 2 map04974 Protein digestion and absorption, 0.00826; map04972 Pancreatic secretion 0.0117
adult gonad sex + morph + food density eggNOG 13 Acid phosphatase A, 0.00452; CUB domain, 0.0139; CwfJ-like, 0.0139; Translation initiation factor IF-2, 0.0139; Mak10, 0.0139; Trypsin-like serine protease, 0.0139; Trypsin, 0.0139; Znf, 0.0139; Allatostatin, 0.0185
juvenile ovaries sex + food density GO 128 GO:0000491(L=7) small nucleolar ribonucleoprotein complex assembly, 0.03146; GO:0000492(L=7) box C/D snoRNP assembly, 0.03146; GO:0001884(L=4) pyrimidine nucleoside binding, 0.03146; GO:0002134(L=4) UTP binding, 0.03146; GO:0002135(L=4) CTP binding, 0.03146; GO:0004360(L=6) glutamine-fructose-6-phosphate transaminase (isomerizing) activity, 0.03146; GO:0007008(L=6) outer mitochondrial membrane organization, 0.03146; GO:0009532(L=6) plastid stroma, 0.03146; GO:0009570(L=7) chloroplast stroma, 0.03146; GO:0009908(L=5) flower development, 0.03146
juvenile ovaries sex + food density KEGG 0
juvenile ovaries sex + food density eggNOG 22 Alpha-2-macroglobulin, 8.29E-03; Alpha-2-macroglobulin receptor, 8.29E-03; Sugar isomerase domain, 8.29E-03; Alpha-macroglobulin thiol-ester catalytic domain, 8.29E-03; DMPK coiled coil domain, 0.0109; Glutamine amidotransferases class-II, 0.0161; Hsp90, 0.0174; CNH domain, 0.0174; DDT DNA-binding domain, 0.0174; p21-Rho-binding domain, 0.0174
adult ovaries sex + morph + food density GO 8 GO:0005319(L=3) lipid transporter activity, 3.24×10^-09; GO:0006801(L=4) superoxide metabolic process, 4.8×10^-09; GO:0006869(L=4) lipid transport, 3.07×10^-08; GO:0072593(L=3) reactive oxygen species metabolic process, 4.60×10^-08; GO:0010876(L=3) lipid localization, 1.33×10^-07; GO:0071702(L=4) organic substance transport, 0.002083; GO:0033036(L=2) macromolecule localization, 0.02463; GO:0006590(L=4) thyroid hormone generation, 0.05054; GO:0006768(L=4) biotin metabolic process, 0.05054; GO:0006797(L=4) polyphosphate metabolic process, 0.05054
adult ovaries sex + morph + food density KEGG 0
adult ovaries sex + morph + food density eggNOG 17 DUF1943, 3.58E-11; Lipoprotein N-terminal Domain, 3.58E-11; Co/Zn superoxide dismutase, 3.58E-11; Vitellogenin, 3.58E-11; von Willebrand factor type D domain, 1.90E-10; MENTAL Cholesterol-capturing domain, 0.00670; DUF3608, 0.00893; Skp1 kinetochore dimerization domain, 0.00893; Skp1 tetramerisation domain, 0.00893; Lipid-binding START domain, 0.0161
juvenile testes sex + food density GO 0
juvenile testes sex + food density KEGG 0
juvenile testes sex + food density eggNOG 0
adult testes sex + morph + food density GO 6 GO:0007278(L=4) pole cell fate determination, 0.03657; GO:0019090(L=5) mitochondrial rRNA export from mitochondrion, 0.03657; GO:0019093(L=4) mitochondrial RNA localization, 0.03657; GO:0030719(L=5) P granule organization, 0.03715; GO:0050790(L=3) regulation of catalytic activity, 0.03715; GO:0065009(L=2) regulation of molecular function, 0.03715; GO:0030695(L=4) GTPase regulator activity, 0.06134; GO:0060589(L=3) nucleoside-triphosphatase regulator activity, 0.06134; GO:0045314(L=9) regulation of compound eye photoreceptor development, 0.0685; GO:0043085(L=4) positive regulation of catalytic activity, 0.0769
adult testes sex + morph + food density KEGG 4 map04360 Axon guidance, 0.0329; map04510 Focal adhesion, 0.0329; map04660 T cell receptor signaling, 0.0329; map04810 Regulation of actin cytoskeleton, 0.0329
adult testes sex + morph + food density eggNOG 9 DUF1630, 0.00711; Myelodysplasia-myeloid leukemia factor 1-interacting protein, 0.00711; P21-Rho-binding domain, 0.00711; Survival motor neuron protein, 0.00711; DENN (AEX-3) domain, 0.00880; dDENN domain, 0.00880; uDENN domain, 0.00880; Tudor domain, 0.0183; RhoGAP, 0.0210
- thorax sex + stage GO 1670 GO:0005622(L=3) intracellular anatomical structure, 1.09×10-104; GO:0071704(L=2) organic substance metabolic process, 4.2×10-97; GO:0043229(L=3) intracellular organelle, 2.40×10-91; GO:0044238(L=2) primary metabolic process, 7.21×10-90; GO:0044237(L=2) cellular metabolic process, 1.21×10-82; GO:0005737(L=4) cytoplasm, 1.87×10^-76; GO:0043227(L=2) membrane-bounded organelle, 8.45×10-66; GO:0043231(L=4) intracellular membrane-bounded organelle, 8.95×10-66; GO:0048856(L=2) anatomical structure development, 1.71×10-57; GO:0007275(L=3) multicellular organism development, 4.9×10-53
- thorax sex + stage KEGG 68 map01100 Metabolic pathways, 5.042×10-44; map05012 Parkinson disease, 3.85×10-14; map01110 Biosynthesis of secondary metabolites, 2.84×10-12; map05010 Alzheimer disease, 2.12×10-8; map00500 Starch and sucrose metabolism, 4.38×10-7; map04145 Phagosome, 5.02×10-7; map00982 Drug metabolism - cytochrome P450, 2.19×10-6; map00980 Metabolism of xenobiotics by cytochrome P450, 3.82×10-6; map00280 Valine, leucine and isoleucine degradation, 7.08×10-6; map03410 Base excision repair, 1.34×10-5
- thorax sex + stage eggNOG 87 TRANS, 1.71×10-120; SIGNAL, 1.21×10-98; COIL, 1.78×10-53; Insect cuticle protein, 9.19×10-21; Major facilitator superfamily, 4.55×10-08; Chitin binding peritrophin-A domain, 3.66×10-07; RNA recognition motif, 3.66×10-07; Chitin-binding domain type 2, 1.08×10-06; Sugar transporter, 2.52×10-06; A1 Propeptide, 5.70×10-06
- ovaries sex + stage GO 9 GO:0006260(L=5) DNA replication, 1.10×10-9</sub; GO:0034061(L=5) DNA polymerase activity, 1.19×10-8</sub; GO:0003964(L=6) RNA-directed DNA polymerase activity, 5.30×10-7</sub; GO:0006278(L=7) RNA-dependent DNA biosynthetic process, 5.30×10-7</sub; GO:0006259(L=4) DNA metabolic process, 1.70×10-6</sub; GO:0016779(L=4) nucleotidyltransferase activity, 0.000765; GO:0015074(L=5) DNA integration, 0.007542; GO:000452(L=8) RNA-DNA hybrid ribonuclease activity, 0.0124; GO:0005319(L=3) lipid transporter activity, 0.0414; GO:0070001(L=5) aspartic-type peptidase activity, 0.0507
- ovaries sex + stage KEGG 0
- ovaries sex + stage eggNOG 5 reverse transcriptase, 1.58×10-8; Rnase H, 0.00235; Endonuclease/Exonuclease/phosphatase family, 0.00240; Integrase, 0.00240; ZnF, 0.0339
- testes sex + stage GO 0
- testes sex + stage KEGG 0
- testes sex + stage eggNOG 0

Return to https://github.com/aphanotus/morphDE